/*
    Copyright 2007-2011 Patrik Jonsson, sunrise@familjenjonsson.org. 
    Based on code by Volker Springel, used with permission.

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Implementation of the multiphase class, which calculates the
    parameters of the Springel & Hernquist 2003 multiphase model.
*/
// $Id$

#include "multiphase.h"
#include <cmath>
#include <iostream>
#include "mcrx-debug.h"

using namespace std;

pair<double, double> 
multiphase::cold_fraction(double density) const
{
  if(density < rho_th)
    return make_pair(0.0,0.0);

  const double tsfr = sqrt(rho_th / density) * t0_star;

  const double A = pow(density / rho_th, -0.8) * A0;

  const double u_h = u_SN / (1 + A) + u_c;

  double ne = 1; // this is just an initial guess and this value seems to work

  const double tcool = c.CoolingTime(u_h, density , ne);

  const double y = tcool>0?
    tsfr / tcool * u_h / (beta * u_SN - (1 - beta) * u_c) : 0.0;

  const double x = y>0? 1 + 1 / (2 * y) - sqrt(1 / y + 1 / (4 * y * y)) : 0.0;

  const double volume_fraction = x*u_c/(u_h - x*u_h + x*u_c);

  assert(x>=0);
  assert(x<=1);
  assert(volume_fraction>=0);
  assert(volume_fraction<=1);

  return make_pair(x,volume_fraction);
}

double multiphase::hotphase_temperature(double density) const
{
  if(density < rho_th)
	return 0.0;
  
  const double factorEVP = pow(density / rho_th, -0.8) * A0;

  const double u_h = u_SN / (1 + factorEVP) + u_c;

  const double hottemp = meanweight_i * (gamma-1) / (k / m_p) * u_h;

  return hottemp;
} 

double multiphase::compute_threshold() const
{
  const double u_h = u_SN / A0;

  const double u4 = 1 / meanweight_i * (1.0 / (gamma-1)) * (k / m_p) * 1.0e4;

  // roughly what it would come out to be with the formula and it's not
  // sensitive to it anyway
  const double dens=2e-23; 
	      
  double ne = 1.0;
  
  const double tcool = c.CoolingTime(u_h, dens, ne);
  
  const double coolrate = u_h / tcool / dens;
  
  const double x = (u_h - u4) / (u_h - u_c);
  
  const double rho_th = x / pow(1 - x, 2) * (beta * u_SN - (1 - beta) * u_c) /
    (t0_star * coolrate);

  return rho_th;
}

